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Abstract 


We  discuss  the  solution  of  sparse  linear  equations  Ky  =  z,  where  £  is 
symmetric  and  indefinite.  Since  exact  solutions  are  not  always  required,  direct 
and  iterative  methods  are  both  of  interest. 1  f 

An  important  direct  method  is  the  Bunch-Parlett  factorization  K  =  U'DU, 
where  U  is  triangular  and  D  is  block-diagonal.  A  sparse  implementation  exists 
in  the  form  of  the  Harwell  code  MA27.  An  appropriate  tttrattvc  method  is  the 
conjugate-gradient-like  algorithm  SYMMLQ,  which  spfves  indefinite  systems 
with  the  aid  of  a  positive-definite  preconditioner. 

For  any  indefinite  matrix  K,  we  show  that  the  VrDU  factorization  can  be 
modified  at  nominal  cost  to  provide  an  “exact”  preconditioner  for  SYMMLQ. 
We  give  code  for  overwriting  the  block-diagonal  matrix  D  produced  by  MA27. 

We  then  study  the  KKT  systems  arising  in  barrier  methods  for  linear  and 
nonlinear  programming,  and  derive  preconditioners  for  use  with  SYMMLQ— 

For  nonlinear  programs  we  suggest  a  preconditioner  based  on  the  “smaller” 
KKT  system  associated  with  variables  that  are  not  near  a  bound.  For  linear 
programs  we  propose  several  preconditioners  based  on  a  square  nonsingular 
matrix  B  that  is  analogous  to  the  basis  matrix  in  the  simplex  method.  The 
aim  is  to  facilitate  solution  of  full  KKT  systems  rather  than  equations  of  the 
form  AulA  Tay  =  r  when  the  latter  become  excessively  ill-conditioned.  /  t  y  _ 
_  J  ' —  oL  °,  Lf.*.  y  ’  i 

Keywords:  indefinite  systems,  preconditioned,  linear  programming,  non¬ 
linear  programming,  numerical  optimization,  barrier  methods,  interior-point 
methods. 
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1.  Introduction 

Symmetric  indefinite  systems  of  linear  equations  arise  in  many  areas  of  scientific 
computation.  We  will  discuss  the  solution  of  sparse  indefinite  systems  Ky  =  z  by 
direct  and  iterative  means. 

The  direct  method  we  have  in  mind  is  the  Bunch-Parlett  factorization  K  = 
UtDU ,  where  U  is  triangular  and  D  is  block-diagonal  with  blocks  of  dimension  1 
or  2  that  may  be  indefinite.  Such  a  factorization  exists  for  any  symmetric  matrix  K 
[BP71].  (We  shall  refer  to  it  as  the  Bunch-Parlett  factorization,  while  noting  that 
the  Bunch-Kaufman  pivoting  strategy  is  preferred  in  practice  [BK77].  The  principal 
sparse  implementation  to  date  is  due  to  Duff  and  Reid  [DR82,DR83]  in  the  Harwell 
code  MA27.  See  also  [DGRST89].) 

The  iterative  method  to  be  discussed  is  the  Paige-Saunders  algorithm  SYMMLQ 
[PS75].  This  is  a  conjugate- gradient-like  method  for  indefinite  systems  that  can 
make  use  of  a  positive-definite  preconditioner. 

1.1.  Preconditioning  indefinite  systems 

One  of  our  aims  is  to  present  a  new  and  simple  result  that  shows  how  to  use  the 
Bunch-Parlett  factorization  of  an  indefinite  matrix  to  construct  an  exact  precondi¬ 
tioner  for  an  iterative  method  such  as  SYMMLQ.  The  intended  use  is  as  follows. 

Given  an  indefinite  system  Ky  =  z  and  a  related  indefinite  matrix  K,  we  expect 
that  the  Bunch-Parlett  factorization  K  =  VTDU  will  be  computed  (or  will  already 
be  available).  We  show  that  D  can  be  changed  cheaply  to  provide  a  positive-definite 
matrix  M  =  UTDU ,  such  that  SYMMLQ  (with  preconditioner  M)  will  solve  Ky  =  z 
in  at  most  two  iterations.  Hence,  M  should  be  a  good  preconditioner  for  the  original 
system  involving  K . 


1.2.  Optimization 

As  a  source  of  indefinite  systems,  we  are  interested  in  barrier  methods  or  interior- 
point  methods  for  solving  linear  and  nonlinear  programs  in  the  following  standard 
form: 

minimize  cTx 

X  (1-1) 

subject  to  Ax  =  b,  l  <  x  <  ti, 

where  A  6  3?mxn,  and 

minimize  F(x) 

1  (1-2) 

subject  to  c(x)  =  0, 


where  F(x)  and  c(x)  have  continuous  first  and  second  derivatives.  We  assume  that 
an  optimal  solution  ( z *,**)  exists,  where  x*  is  a  set  of  Lagrange  multipliers  for  the 
constraints  Ax  =  b  or  c(z)  =  0. 
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1.3.  KKT  systems 

When  barrier  or  interior-point  methods  are  applied  to  these  optimization  problems, 
the  Karush- Kuhn- Tucker  optimality  conditions  lead  to  a  set  of  equations  of  the  form 

*■(:')■  «•» 

whose  solution  usually  dominates  the  total  computation.  The  vectors  Ax  and  Ax 
are  used  to  update  the  estimates  of  x*  and  x* . 

For  quadratic  programs  or  general  nonlinear  programs,  H  is  typically  a  general 
sparse  matrix  like  A,  and  it  is  natural  to  solve  the  KKT  system  as  it  stands.  The 
Harwell  code  MA27  has  been  used  in  this  context  by  several  authors,  including  Gill  et 
at.  [GMSTW86]  and  Turner  [Tur87,Tur90]  for  sparse  linear  programs,  by  Ponceleon 
[Pon90]  for  sparse  linear  and  quadratic  programs,  and  by  Burchett  [Bur88]  for  some 
large  nonlinear  programs  arising  in  the  electric  power  industry. 


1.4.  Avoiding  AD2AT 

If  H  is  known  to  be  nonsingular,  it  is  common  practice  to  use  it  as  a  block  pivot 
and  solve  (1.3)  according  to  the  range-space  equations  of  optimization: 

AH-1y4TAx  =  AH~xg  +  r,  H Ax  =  ATAx  -  g. 

For  linear  programs  this  is  particularly  attractive,  since  H  is  then  a  positive  diagonal 
matrix.  For  example,  in  a  typical  primal  barrier  method,  H  —  pD~2  where  D  is 
diagonal  and  p  is  the  barrier  parameter  {p  >  0)  [GMSTW86].  The  range-space 
equations  reduce  to 

AD2ATAir  =  AD2g  +  pr ,  Ax  ■=  ~D2(ATArc  -  g),  (1.4) 

and  most  of  the  work  lies  in  solving  the  system  involving  AD2AT.  When  r  =  0,  the 
numerical  properties  may  be  improved  by  noting  that  the  equation  for  Ax  reduces 
to  the  least-squares  problem 

min  \\Dg  -  DA7 Ax||2.  (1.5) 

A»r 

However,  it  is  important  to  observe  that  the  range-space  equations  map  not  give  a 
stable  method  for  solving  the  KKT  system  if  H  is  ill-conditioned. 


1.5.  Example 
Let 


fl  1 

\ 

A  = 

1 

1 

,  H  = 

l  1 

L 

\ 

f  v/£  \ 

,  X  = 

1 

l  1 
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where  /*  <  1,  and  consider  the  KKT  system  (1.3).  This  would  arise  when  a  primal 
barrier  method  is  applied  to  a  3  x  4  LP  problem  (1.1)  having  /  =  0,  u  =  oo,  when  x  is 
the  current  estimate  of  x*  and  /*  is  the  current  barrier  parameter.  Thus  H  =  pD~2, 
where  D  =  diag(xj). 

The  condition  numbers  of  interest  are  cond(A')  «  6  (independent  of  p)  and 
cond(AH~x At)  =  cond(AD2AT)  «  0.5//*. 1  The  latter  becomes  increasingly  large 
as  a  solution  is  approached  (/*  — ♦  0),  even  though  K  and  the  original  linear  program 
are  very  well-conditioned. 

Similar  examples  are  easily  constructed.  (Indeed,  K  can  be  well-conditioned 
even  if  H  is  singular.)  Thus,  we  advocate  direct  or  iterative  solution  of  the  full 
KKT  system  (1.3)  even  for  linear  programs,  rather  than  (1.4)  or  (1.5)  according  to 
current  practice. 

Gay  [Gay89,  pp.  16-17]  has  already  drawn  attention  to  the  lurking  numerical 
difficulties  and  suggests  a  middle  ground  of  working  with  AD2AT  as  long  as  possible, 
then  switching  to  a  more  robust  alternative  such  as  direct  solves  with  K. 

1.6.  Iterative  methods  and  preconditioning 

The  KKT  systems  we  are  concerned  with  arise  when  Newton’s  method  is  applied 
to  the  nonlinear  equations  defining  optimality  conditions  for  barrier  subproblems; 
see  Section  3.  In  this  context,  there  are  not  many  KKT  systems  to  be  solved 
(compared  to  those  in  active-set  methods),  the  systems  need  not  be  solved  exactly 
(Dembo,  Eisenstat  and  Steihaug  [DES82]),  and  the  KKT  matrix  eventually  does 
not  change  significantly.  It  is  therefore  appropriate  to  consider  iterative  methods 
and  preconditioners  for  the  indefinite  matrix  K. 

Previous  work  on  preconditioning  for  interior-point  methods  has  focussed  on  the 
LP  case  and  the  Schur-complement  matrix  AD2AT.  Most  approaches  have  used 
approximate  Cholesky  factors  of  AD2AT\  for  example,  [GMSTW86,Kar87,KR88, 
Meh89a].  Exact  LU  factors  of  DAT  have  also  been  investigated  [GMS89]. 

The  success  of  preconditioned  conjugate-gradient  methods  in  this  context  lends 
added  promise  to  our  proposed  use  of  the  much  better-conditioned  KKT  systems, 
now  that  it  is  known  how  to  precondition  indefinite  systems. 

1.7.  Summary 

In  Section  2  we  consider  general  indefinite  systems  and  derive  a  preconditioner 
from  the  Bunch-Parlett  factorization.  In  Section  3  we  consider  barrier  methods  for 
nonlinear  programs,  and  propose  factorizing  just  part  of  the  KKT  system  to  obtain 
a  preconditioner  for  the  whole  system. 

Sections  4-6  deal  with  the  LP  case.  In  Section  4  we  propose  three  preconditioners 
based  on  LU  factors  of  a  square  nonsingular  matrix  B  (analogous  to  the  basis  in 
the  simplex  method).  Section  5  discusses  some  practical  difficulties.  Section  6  gives 
numerical  results  on  the  condition  numbers  of  K  and  A  D2AT  in  a  typical  sequence 
of  barrier  subproblems,  and  compares  the  preconditioned  systems  C~XKC~T  for 
several  preconditioners  CCT. 

'We  use  the  spectral  condition  number,  cond(A')  =  ||A'-,||j||A'||j. 
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2.  Preconditioning  Indefinite  Systems 

Let  K  be  any  symmetric  nonsingular  matrix,  and  let  M  be  a  given  positive-definite 
matrix.  Also,  let  “products  with  K"  mean  matrix-vector  products  of  the  form 
u  =  Kv,  and  “solves  with  Af”  mean  solution  of  linear  systems  of  the  form  Mx  =  y. 

The  Paige-Saunders  algorithm  as  implemented  in  SYMMLQ  [PS75]  may  be  used 
to  solve  Ky  =  z  even  if  K  is  indefinite.  As  with  other  conjugate-gradient-like 
algorithms,  the  matrix  is  represented  by  a  procedure  for  computing  products  with 
K  (those  generated  by  the  symmetric  Lanczos  process). 

The  first  steps  towards  accelerating  the  conve’gence  of  this  algorithm  were  taken 
by  Szyld  and  Widlund  [SW78].  Given  a  positive-definite  matrix  Af  as  precondi¬ 
tioner,  their  algorithm  used  solves  with  Af  in  the  normal  way,  but  was  unconven¬ 
tional  in  also  requiring  products  with  M. 

Subsequently,  a  variant  of  SYMMLQ  was  developed  that  requires  only  solves  with 
Af  [Sau79].  To  solve  Ky  =  z,  this  variant  regards  the  preconditioner  as  having  the 
form  Af  =  CCTand  implicitly  applies  the  Paige-Saunders  algorithm  to  the  system 

C~lKC~Tw  =  C_12, 

accumulating  approximations  to  the  solution  y  =  C~Tw  (without  approximating  w, 
which  isn’t  needed).  An  implementation  is  available  from  the  misc  chapter  of  netlib 
[DG85]. 

2.1.  Use  of  the  Bunch-Parlett  factorization 

Given  any  symmetric  nonsingular  matrix  K ,  there  exists  a  factorization  of  the  form 

K  =  PtUtDUP, 

where  P  is  a  permutation,  U  is  upper  triangular,  and  D  is  block-diagonal  with 
blocks  of  dimension  1  or  2  [BP71j.  If  K  is  indefinite,  some  of  the  blocks  of  D  will 
have  negative  eigenvalues.  Let  the  eigensystem  of  D  be 

D  =  QAQt,  A  =  diag(Aj), 

and  let 

D  =  QAQt,  A  =  diag(|Aj|), 

be  a  closely  related  positive-definite  matrix  that  can  be  obtained  at  minimal  cost. 
If  we  define  C  =  PTUTD 1^2,  it  is  easily  verified  that 

K  =  C~'KC-t  =  diag(A;/|A,|)  =  diag(±l). 

This  means  that  the  “perfect”  preconditioner  for  I\  is  the  matrix 

M  =  CCT  =  PTUTDUP , 

since  the  “preconditioned”  matrix  Ii  has  at  most  two  distinct  eigenvalues  and  the 
Paige-Saunders  algorithm  converges  in  at  most  two  iterations. 

In  practice,  M  will  be  computed  from  the  Bunch-Parlett  factorization  of  an 
approximation  to  K. 
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2.2.  Modification  of  D  from  MA27 

The  block-diagonal  matrix  D  is  packed  in  the  MA27  data  structure  as  a  sequence  of 
matrices  of  the  form 

(“)  “d  (£7)' 

In  the  lxl  case,  we  do  nothing  if  a  >  0;  otherwise  we  reverse  its  sign.  In  the  2x2 
case,  we  do  nothing  if  07  >  (32\  otherwise  we  compute  the  eigensystem  in  the  form 

(; >,)(:-:)■ 

where  <?  +  s2  =  1.  We  then  form  the  positive-definite  matrix 


and  overwrite  the  appropriate  three  locations  of  MA27’s  storage. 

The  techniques  for  computing  the  2x2  eigensystem  are  central  to  Jacobi’s 
method  for  the  symmetric  eigenvalue  problem.  They  were  developed  by  Rutishauser 
[Rut66].  We  have  followed  the  description  in  Golub  and  Van  Loan  [GV89,  p.  446], 
with  minor  changes  to  work  with  symmetric  plane  rotations. 

A  subroutine  for  modifying  the  D  computed  by  MA27  is  given  in  the  Appendix. 

2.3.  Aasen’s  method 

In  general,  Aasen’s  tridiagonalization  method  [Aas71]  is  considered  competitive  with 
the  Bunch-Kaufman  approach  [BK77]  for  solving  dense  indefinite  systems.  Aasen’s 
method  computes  a  factorization  of  the  form  K  —  UtTU  where  T  is  tridiagonal. 

We  do  not  know  of  a  sparse  implementation,  but  in  any  event  we  note  that  it 
would  not  be  ideal  for  producing  a  preconditioner  in  the  manner  described  above, 
since  the  eigensystem  for  T  would  involve  far  more  work  than  for  the  block-diagonal 
D  of  the  Bunch-Parlett  factorization. 

On  the  other  hand,  we  could  compute  a  (very  special)  Bunch-Parlett  factoriza¬ 
tion  of  T  and  modify  the  associated  D  as  described  above. 
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3.  Barrier  Subproblems 

We  return  now  to  the  optimization  problems  (1.1)— (1.2).  In  barrier  methods,  the 
bounds  l  <  x  <  u  are  absorbed  into  the  objective  function  and  we  solve  a  sequence 
of  perturbed  subproblems,  typically  of  the  form 


minimize  F„(x)  =  F(x )  -  /i£(ln (ij  -  lj)  +  ln(uj  -  xj)) 

j= l  (3.1) 

subject  to  c(x)  =  0, 

where  the  barrier  parameter  /z  takes  decreasing  positive  values  that  are  eventually 
very  small.  If  lj  =  -oo  or  Uj  =  oo  for  some  j,  the  corresponding  terms  ln(xj  -  lj) 
or  ln(«j  —  Xj)  are  omitted.  If  Xj  has  no  bounds,  both  terms  may  be  omitted. 

The  following  quantities  are  needed: 


Lm(x,x) 

=  F„(x)  -  7TTc(x) 

9n(x) 

=  VF„(x) 

9l(x,tt) 

=  5M(x)  -  A(x)Tv 

Hl(x,x) 

=  V2/;(x)-S>,V2cf(x) 

A(x) 

=  Vc(x) 

the  Lagrangian  function, 
the  gradient  of  the  barrier  function, 
the  gradient  of  the  Lagrangian, 
the  Hessian  of  the  Lagrangian,  and 
the  Jacobian  of  the  constraints. 


For  convenience  we  assume  that  A{x)  £  J?mXn  has  full  row  rank  m,  and  that  the 
scaling  of  the  problem  is  reasonable,  so  that  ||j4(x)||  «  1. 


3.1.  Newton’s  method  and  the  KKT  system 

The  optimality  conditions  for  (3.1)  are  the  nonlinear  equations 

gdx,  tt)  =  0,  (3.2) 

c(x)  =  0.  (3.3) 

Newton’s  method  may  be  applied  directly,  or  to  some  equivalent  system.  Given 
suitable  initial  values  for  the  primal  and  dual  variables  (x,tt),  the  key  set  of  equations 
for  generating  a  search  direction  is  the  KKT  system 

('■  ")(-£)-(:)■ 

where  the  KKT  matrix  and  right-hand  side  are  evaluated  at  the  current  point  (x,tt). 
A  positive  steplength  a  is  then  chosen  to  reduce  some  measure  of  the  size  of  the 
right-hand  side  (<ft.,c),  and  the  variables  are  updated  according  to  x  «—  x  +  a  Ax, 
x  «—  jr  4-  qAtt.  (Sometimes  a  different  a  may  be  used  for  x  and  7r.) 
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3.2.  A  preconditioner  for  K 

In  general,  some  of  the  variables  converge  to  values  near  their  upper  or  lower  bounds. 
For  such  variables  Xj,  the  Hessian  HL  includes  on  its  diagonal  a  term  that  becomes 
very  large:  n/(xj  -  lj)2  or  n/(uj  -  Xj)2,  which  are  0(1/ fi).  Let  the  KKT  matrix  be 
partitioned  accordingly: 


where  K\  is  the  part  of  HL  associated  with  variables  near  a  bound,  and  Ki  looks  like 
a  smaller  KKT  system  associated  with  the  remaining  variables.  This  partitioning 
is  crucial  to  the  sensitivity  analysis  in  [Pon90].  Of  course,  the  partition  depends 
on  the  measure  of  closeness  to  a  bound,  but  it  is  not  critical  here  except  that  the 
dimension  of  Ii\  should  not  exceed  n  -  m. 

One  possible  approximation  to  A*  is 


where  D i  is  a  diagonal  matrix  containing  the  diagonals  of  A'i,  which  by  construction 
are  large  and  positive.  Applying  the  method  of  Section  2,  we  can  now  obtain  a 
positive-definite  preconditioner  for  A'  as  follows: 

M=(D'  (3.7| 

where  D2  is  obtained  from  Z?2  at  nominal  cost. 

3.3.  Discussion 

In  broad  terms,  we  need  to  estimate  which  variables  are  going  to  be  “free”  (away 
from  their  bounds)  at  a  solution.  If  m  <  n,  the  KKT  system  A'2  associated  with 
the  free  variables  may  be  much  smaller  than  the  whole  of  K,  and  the  cost  of  the 
Bunch-Parlett  factorization  of  A'2  may  be  acceptably  low. 

For  the  early  iterations  of  Newton’s  method,  the  estimate  of  A'2  will  usually  be 
poor,  and  the  diagonal  term  D\  will  not  be  particularly  large.  However,  following  the 
inexact-Newton  approach  [DES82],  only  approximate  solutions  to  the  KKT  system 
are  needed,  and  the  iterative  solver  need  not  perform  many  iterations. 

As  the  Newton  iterations  converge  and  the  partition  (3.5)  becomes  more  sharply 
defined,  the  preconditioner  should  become  increasingly  powerful  and  produce  the 
increasingly  accurate  solutions  required  at  an  acceptable  cost. 
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4.  Preconditioners  for  Linear  Programming 

For  linear  programs  the  structure  of  the  partitioned  KKT  system  (3.5)  can  be  inves¬ 
tigated  more  closely,  given  that  optimal  solutions  are  often  associated  with  a  vertex 
of  the  feasible  region.  We  partition  the  constraint  matrix  into  the  form  A  =  (N  B  ), 
where  B  is  square  and  nonsingular,  and  N  in  some  sense  corresponds  to  the  n  —  m 
variables  that  are  closest  to  a  bound. 

The  Hessian  for  the  barrier  function  is  a  diagonal  matrix  H,  which  we  partition 
as  H  =  diag (HS,HB).  The  KKT  system  is  then 

(  Hn  Nt  \ 

K  =  Hb  Bt  . 

\  N  B 

As  convergence  occurs,  the  diagonals  of  HN  — ►  oc  (and  in  general  cond(A')  — »  oo). 
In  degenerate  cases,  some  diagonals  of  HB  may  also  become  very  large. 

In  various  primal,  dual  and  primal-dual  interior-point  algorithms  for  LP,  sim¬ 
ilar  matrices  K  arise  with  varying  definitions  of  H  (e.g.,  [Meg86,KMY88,LMS89, 
Meh89a,Meh90]).  The  discussion  hereafter  applies  to  all  such  methods. 

In  the  following  sections  we  introduce  a  series  of  preconditioners  of  the  form 
M  =  CCT.  To  improve  the  convergence  of  SYMMLQ,  the  transformed  matrices  K  = 
C~*KC~T  should  have  a  better  condition  than  K  or  a  more  favorable  distribution 
of  eigenvalues  (clustered  near  ±1).  We  make  use  of  the  quantities 

V  =  B~tHbB-\  W  =  NH^l/2, 

and  are  motivated  by  the  fact  that  V  — ■  0  and  W  — *■  0  in  nondegenerate  cases.  The 
effects  of  degeneracy  are  discussed  later. 


4.1.  The  preconditioner  A/j 


The  first  preconditioner  is  diagonal  and  is  intended  to  eliminate  the  large  diagonals 
of  K: 


(  hn 


Mx  =  CXC{  = 


V 


\ 


/ 


/  I 


/?!  =  cr'/vcf7  = 


WT  \ 


hb  Bt 
\W  B 


(4.1) 


(4.2) 


With  diagonal  preconditioning,  there  is  no  loss  of  precision  in  recovering  solutions 
for  the  original  system.  Thus  as  HN  becomes  large,  the  preconditioned  matrix  h\ 
tends  to  represent  the  true  sensitivity  of  the  KKT  system  with  regard  to  solving 
linear  equations. 

We  will  use  Kx  later  for  comparing  condition  numbers 
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4.2.  The  preconditioner  M2 

The  second  preconditioner  is  block  diagonal: 

Hn 

m2  =  C2C2  = 


btb 


( I 


A',  =  C;xKC;t  = 


I 

f 

WT 

/ 


(4.3) 


(4.4) 


V 

\  W  I 

Since  V  and  IV  tend  to  become  small,  M2  tends  towards  being  an  exact  precon¬ 
ditioner  for  K.  We  see  that  a  Bunch-Parlett  factorization  is  no  longer  needed.  In 
order  to  solve  systems  involving  M2 ,  we  may  use  any  sparse  factorization  of  B  or 
BT. 


4.3.  The  preconditioner  M3 

The  third  preconditioner  is  designed  to  eliminate  the  submatrix  V  in  (4.4),  for 
degenerate  cases  where  V  is  not  adequately  small: 


Mj  —  C3C^,  Co  — 


H'J2 


Br  \HbB~x  , 


(4.5) 


K,  =  C;'KCVt  = 


(  / 
-\vw 


.lw^v  WT 


(4.6) 

\  W  I 

The  off-diagonal  term  in  (4.5)  can  be  derived  by  observing  that  for  a  KKT  matrix 
of  the  form 


-(:")■ 


B  square. 


we  would  like  M  =  CCT  to  satisfy 

C"''C‘T  =  (, 

or  equivalently,  CJCT  =  K.  Letting  C  be  of  the  form 


we  find  that  E  should  satisfy  EB  +  BTET  =  H.  The  simplest  choice  is  then  to  set 
E  =  \HB~l. 

Though  V  has  been  eliminated,  we  have  now  introduced  the  term  -%VW,  and 
solves  with  M3  cost  twice  as  much  as  solves  with  M2.  The  expected  benefit  is  that 
^VW  should  be  smaller  than  V  itself. 
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4.4.  The  preconditioner  M4 

The  fourth  preconditioner  also  eliminates  V ,  using  the  factorization  BT  =  LU ,  where 
we  intend  that  L  be  well-conditioned: 

( H'J2 

M,  =  C\Cj,  C\  =  I  L  ,  (4.7) 

V  uT  ) 

f  I  -\WTV  WT  N 

I(4  =  C41KC4T  =  -\VW  I  .  (4.8) 

\  W  I 

where 

v  =  l-xhbl~t,  W  =  U~tNHsXI2. 

As  before,  letting  C  be  of  the  form 


and  requiring  CJCr—  K,  we  find  that  ELT+  LET  =  H ,  and  we  take  E  =  %HL~T. 

Solves  with  M4  are  cheaper  than  with  M3.  Comparing  (4.6)  and  (4.8),  a  further 
advantage  is  that  VW  -  UVW  tends  to  be  smaller  than  VW,  although  W  —  U~TW 
is  probably  larger  than  W. 

4.5.  The  preconditioner  MD 

We  mention  one  further  diagonal  preconditioner  that  has  appeared  implicitly  in  the 
literature  for  the  case  H  =  (iD~2  with  D  diagonal.  It  does  not  depend  on  the  N-B 
partitioning,  and  gives  a  transformed  system  that  does  not  involve  /r. 

MD  =  CDCj=  ^ D~ 2  ^_XJ  (4.9) 

■  c°',<c°T -  ( L  DAT ) •  ,4io) 

The  matrix  IiD  is  associated  with  weighted  least-squares  problems  of  the  form  (1.5), 
as  discussed  in  [GMSTW86].  Turner  [Tur87,Tur90]  has  investigated  the  use  of  MA27 
to  obtain  exact  factors  of  both  K  and  I\D.  An  important  practical  observation  was 
that  MA27  produced  much  sparser  factors  for  1\ D  than  for  1\. 

Unfortunately,  the  numerical  examples  in  Section  6  show  that  KD  has  essentially 
the  same  condition  as  AD2Ar,  which  tends  to  be  much  more  ill-conditioned  than 
Kx  (4.2)  .  We  therefore  cannot  recommend  the  use  of  KD. 

For  the  sake  of  both  direct  and  iterative  methods  for  solving  KKT  systems,  it  is 
hoped  that  further  development  of  MA27  will  result  in  greatly  improved  sparsity  in 
the  factors  of  I\  and/or  A',.  At  the  time  of  writing,  a  new  code  MA47  holds  much 
promise  (see  Duff  et  al.  [DGRST89]). 
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4.6.  Regularizing  K  and  AD2AT 

Since  A  often  does  not  have  full  row  rank,  it  is  important  to  include  a  regularization 
parameter  6  >  0  in  the  KKT  system.  Thus  (1.3)  becomes 

C 

Systems  of  this  type  have  been  studied  in  the  context  of  sequential  quadratic  pro¬ 
gramming  by  Murray  [Mur69],  Biggs  [Big75]  and  Gould  [Gou86]. 

In  practice,  a  wide  range  of  values  of  6  may  be  used  without  inhibiting  conver¬ 
gence,  particularly  with  methods  that  do  not  maintain  primal  feasibility  (AAx  =  0). 
For  example,  we  would  recommend  values  in  the  range  10-8  <  6  <  10-4  on  a  ma¬ 
chine  with  about  16  digits  of  precision,  assuming  ||.A||  ss  1. 

Note  that  the  corresponding  system  (1.4)  becomes 

(AD2AT+  ft6I)Air  =  AD2g  +  fir.  (4.12) 

When  /i  is  as  small  as  10-10  (say),  one  would  have  to  choose  rather  a  large  6  (say 
6  >  10-2)  to  achieve  any  degree  of  regularization  of  AD2AT.  This  constitutes  a  large 
perturbation  to  the  underlying  KKT  system  (4.11). 

In  other  words,  a  much  smaller  6  is  sufficient  to  regularize  (4.11)  than  (4.12). 
Thus,  KKT  systems  again  show  an  advantage  over  AD2AT. 

With  regard  to  the  preconditioners,  6  introduces  terms  —61,  —61,  -6U~TU~l 
into  the  bottom  corner  of  Ii'2 ,  I(3,  K4  respectively.  For  Ii4  it  appears  that  6  must 
be  chosen  quite  small  and  that  the  choice  of  B  must  be  flexible  enough  to  prevent 
U  from  being  excessively  ill-conditioned  (see  Section  5.3). 


IS 


5.  Use  of  LU  Factors 

For  linear  programs,  the  “small”  KKT  matrix  in  (3.5)  is  of  the  form 


As  in  the  general  nonlinear  case  we  could  obtain  a  preconditioner  from  a  Bunch- 
Parlett  factorization  of  K2,  and  in  practice  this  may  prove  to  be  a  good  approach. 

The  preconditioners  Af2,  M3  and  M4  were  derived  on  the  assumption  that  it 
should  be  cheaper  to  compute  sparse  factors  of  just  the  matrix  B.  We  propose  to 
use  the  package  LUSOL  [GMSW87]  to  obtain  BT  =  LU ,  where  L  is  a  permuted 
lower  triangle  with  unit  diagonals.  A  user-defined  tolerance  limits  the  size  of  the 
off-diagonals  of  L  (typically  to  5,  10  or  100),  thereby  limiting  the  condition  of  L  as 
required. 

5,1.  Choice  of  B 

One  of  the  main  practical  difficulties  will  be  in  choosing  a  “good”  square  matrix  B  at 
each  stage.  The  current  values  of  x  and /or  the  estimated  reduced  costs  z  —  c  —  ATir 
should  provide  some  guidance.  For  example,  the  diagonal  matrix  H  is  defined  in 
terms  of  these  quantities,  and  the  smallest  m  +  s  diagonals  of  H  could  be  used  to 
pinpoint  a  submatrix  A  of  A  (for  some  moderate  $  >  0).  LUSOL  could  then  be  used 
to  obtain  a  rectangular  factorization  AT=  LU.  The  first  m  pivot  rows  and  columns 
may  suggest  a  suitable  B. 

Alternative  approaches  to  choosing  B  have  been  suggested  by  Gay  [Gay89],  Tapia 
and  Zhang  [TZ89],  Mehrotra  [Meh89b]  and  others.  These  remain  to  be  explored. 


5.2.  The  effects  of  degeneracy  on  V  and  W 

In  general,  primal  degeneracy  will  mean  that  certain  elements  of  HB  do  not  tend  to 
zero,  so  that  not  all  of  V  or  V  will  become  small.  Similarly,  dual  degeneracy  will 
mean  that  certain  elements  of  HN  will  not  become  large,  and  not  all  of  W  or  W  will 
become  small. 

The  main  effect  is  that  the  preconditioners  will  be  less  “exact”.  Either  form  of 
degeneracy  is  likely  to  increase  the  number  of  SYMMLQ  iterations  required. 


5.3.  Singular  systems 

Whatever  the  method  for  choosing  a  square  B ,  it  is  probable  that  B  will  be  singular 
(since  in  many  practical  cases,  A  does  not  have  full  row  rank).  At  present  we  propose 
to  rely  on  the  fact  that  LUSOL  will  compute  a  stable  singular  factorization  of  the 
form 


fiT=  ( i! )  ( 1,1  v *)• 


u 
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and  the  solve  procedures  will  treat  this  as  if  it  were  the  factorization  of  a  nonsingular 
matrix 

"-U ,)("  ’;)■ 

User-defined  tolerances  determine  how  ill-conditioned  U\  is  allowed  to  be  (and  hence 
determine  its  dimension). 

Alternatively,  we  may  use  the  factorization  B7  =  L\U\  to  transform  most  of  K 
as  already  described.  Certain  rows  of  A  will  not  be  transformed  in  the  preferred  way, 
and  again  the  effect  will  be  to  increase  the  number  of  SYMMLQ  iterations  required. 

6.  Numerical  Examples  for  the  LP  Case 

Here  we  investigate  the  effect  of  the  preconditioners  described  in  Section  4.  For  test 
purposes  we  have  used  MATLAB  [MLB87]  to  implement  a  primal-dual  interior- 
point  algorithm  for  the  standard  LP  problem  min  cTx  subject  to  Ax  =  b,  x  >  0. 
The  linear  system  to  be  solved  each  iteration  is 

(  A  -SI  )  (-11  M"r)'  (6  ‘> 

where  H  =  X~1Z,  X  =  diag(ij),  Z  =  diag (zj),  r  —  b  —  Ax,  g  =  c  -  A7 ir  -  pX-1e, 
and  e  is  the  vector  of  ones.  The  search  direction  for  z  is  Az  =  X~l(fie  -  ZAx)  -  z. 

The  rows  and  columns  of  A  were  scaled  to  give  ||A||  1.  The  starting  values 

were  x  =  e,  z  —  e,  ir  =  0  (so  that  H  =  I  initially),  and  6  was  fixed  at  10~8,  with  the 
machine  precision  on  a  DEC  VAX  system  being  around  16  digits.  The  parameter  n 
was  reduced  every  iteration  according  to  the  steplengths  for  x  and  z:  p  <—  p  -  a^p, 
where  =  min(a*,a2,0.99)  and  ar,  az  were  limited  in  the  usual  way  to  be  at 
most  1  or  0.99  times  the  step  to  the  boundaries  x  >  0,  z  >  0  respectively.  See 
[KMY88,MMS89,LMS89,Meh90]  for  related  details. 

Condition  numbers  of  various  matrices  were  obtained  using  MATLAB’s  function 
rcond.  The  square  matrices  B  for  the  preconditioners  of  Section  4  were  obtained 
from  the  columns  of  A  for  which  Hjj  <  20.  The  diagonals  Hn  were  first  sorted  and 
up  to  1.2m  of  the  smallest  were  used  to  select  a  rectangular  matrix  A  from  A.  In 
practice,  a  sparse  LU  factorization  of  A  or  AT  would  extract  a  full-rank  submatrix, 
but  here  we  used  MATLAB’s  function  qr(A )  to  elicit  a  full- rank  set  of  columns  (via  a 
QR  factorization  with  column  interchanges),  and  a  second  QR  factorization  of  part 
of  ATto  pinpoint  a  full-rank  set  of  rows.  The  dimension  of  the  resulting  matrix  B  is 
generally  less  than  m.  The  “rank”  was  determined  from  the  first  QR  factorization 
by  requiring  the  diagonals  of  R  to  be  greater  than  10-4. 

6.1.  A  nondegenerate  example 

To  illustrate  ideal  behavior  of  the  preconditioners,  we  chose  a  nondegenerate  problem 
ezpl  [Bla82]  in  which  A  is  10  by  17  (including  10  unit  columns  associated  with  slack 
variables).  The  lack  of  primal  or  dual  degeneracy  means  that  near  a  solution,  m  =  10 
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diagonals  of  H  are  substantially  less  than  1  and  n-m  =  7  diagonals  are  significantly 
greater  than  1.  The  choice  of  B  is  ultimately  clear-cut. 

Table  1  lists  various  condition  numbers  for  each  iteration  of  the  primal-dual 
algorithm.  For  interest  we  include  AD2AT  and  KD,  which  were  defined  in  terms 
of  D  —  X  =  diag(ij)  (see  Section  4.5)  and  incorporated  the  same  regularization  6 
(Section  4.6).  It  may  be  seen  that  both  AD2AT+fi6I  and  K D  become  increasingly  ill- 
conditioned  in  step  with  K ,  in  contrast  to  the  “meaningful”  condition  of  K  reflected 
by  Kl  (in  which  the  large  diagonals  of  H  have  been  scaled  to  1). 

The  preconditioned  systems  K2 ,  K 3  and  K4  show  an  increasing  though  appar¬ 
ently  mild  improvement  over  Kl.  Their  effectiveness  depends  on  the  choice  of  B 
and  whether  or  not  it  has  dimension  m.  The  column  labeled  aB  rank-def”  records 
the  corresponding  rank  deficiency.  The  conditions  of  B,  L  and  U  were  less  than  25, 
7  and  40  respectively  for  all  iterations. 

Low  conditions  are  always  a  good  sign,  but  high  ones  tell  an  incomplete  story. 
Figure  1  shows  more  clearly  the  increasing  improvement  of  the  preconditioners  M2, 
M3,  M4  in  terms  of  the  clustering  of  the  eigenvalues  of  I(2,  K 3,  K4  around  ±1. 
The  KKT  systems  have  dimension  m  +  n  =  27.  Eigenvalues  in  the  range  (—5,5)  are 
plotted  exactly;  the  remainder  are  compressed  into  the  ranges  (-6,-5)  and  (5,6). 
Thus,  Ii2  has  one  or  two  eigenvalues  greater  than  5  for  the  first  eight  iterations, 
whereas  Ii3  has  its  eigenvalues  inside  (-5,5)  at  all  times.  (The  vertical  axis  is 
“iteration  number”  shifted  by  1  for  /v2,  14  for  /v3,  and  27  for  K4.  Each  horizontal 
line  gives  the  spectrum  of  one  of  these  matrices  at  the  corresponding  iteration.) 

It  is  evident  from  Figure  1  that  I(3  and  K4  have  more  favorable  eigenvalue 
distributions  than  I\2,  and  that  1\4  is  marginally  better  than  A'3,  the  main  benefit 
being  that  it  is  more  cheaply  obtained.  There  is  a  striking  absence  of  eigenvalues  in 
the  range  (- 1  +  /?,  -/3)  for  some  small  /?,  though  we  have  no  immediate  explanation. 
This  range  broadens  to  (  —  1  +  /?,  1  -  fl)  for  all  systems  at  the  final  iteration,  as  we 
may  expect. 

6.2.  A  more  typical  example 

Table  2  and  Figure  2  give  similar  results  for  the  well  known  problem  afiro  [Gay85]. 
The  matrix  A  is  27  by  51,  including  19  slack  columns.  We  see  that  AD2AT  +  pSI 
and  Kd  again  become  extremely  ill-conditioned  in  step  with  K. 

The  KKT  systems  have  dimension  78.  As  before  there  is  a  clear  division  between 
large  and  small  diagonals  of  H  near  a  solution,  but  in  this  case  only  m  -  5  are 
substantially  smaller  than  one.  The  rank  of  the  corresponding  columns  of  A  is 
m  -  7,  consistent  with  B' s  final  rank  deficiency  of  7.  The  conditions  of  B,  L  and  U 
were  again  low:  less  than  35,  13  and  34  respectively. 

It  is  encouraging  to  observe  that  Figure  2  is  qualitatively  similar  to  Figure  1 
in  spite  of  the  rank  deficiency  in  B.  The  main  difference  is  two  eigenvalues  close 
to  zero  on  the  last  iteration,  in  keeping  with  the  difference  between  m  -  5  and 
m  -  7.  Similarly  for  the  second  to  last  iteration.  We  can  expect  a  low  number  of 
SYMMLQ  iterations  will  be  required  as  the  barrier  algorithm  converges,  as  in  the 
ideal  nondegenerate  case. 
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Jfc 

AD2At 

kd 

K 

Ki 

k2 

*3 

*4 

B  rank-def 

1 

1.9e-2 

1.2el 

l.lel 

l.lel 

l.lel 

3. lei 

2.6el 

l.lel 

2 

4.5e-3 

6.5el 

3. lei 

6.8e2 

1.9el 

5.5el 

2.7el 

2.7el 

3 

4.8e-4 

1.2e3 

2.2e3 

4.4e4 

6.9el 

4.6e2 

5.7el 

1.6el 

4 

1.3e-4 

1.7e5 

3.1e5 

1.0e6 

4.2e3 

4.4e3 

1.5e3 

1.6e3 

2 

5 

3.2e-5 

2.9e5 

5.3e5 

1.3e6 

7.7e2 

2.3e3 

5.4e2 

5.6e2 

1 

6 

1.5e-5 

4.0e5 

5.8e5 

3.1e5 

6.3e2 

2.8e3 

7.4e2 

7.5e2 

1 

7 

4.6e-6 

4.0e5 

5.8e5 

1.6e5 

1.2e2 

3.0e2 

l.le2 

l.le2 

1 

8 

1.6e-6 

4.7e5 

6.3e5 

3.0e5 

3.9el 

1.4e2 

1.9el 

8.4e0 

9 

1.4e-7 

9.1e5 

l.le6 

4.9e5 

2.0el 

4.7e0 

1.3e0 

1.5e0 

10 

7.8e-9 

8.4e5 

1.0e6 

3.5e6 

1.7el 

l.leO 

l.leO 

1.3e0 

Table  1:  Condition  numbers  for  problem  exp  1 
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* 

AD2At 

Rd 

K 

Rx 

Ri 

R  3 

R  A 

B  rank-def 

1 

2.6e-2 

6.0el 

2.2el 

3.3el 

2.3el 

1.3e2 

1.8e2 

9.8el 

1 

2 

9.9e-3 

3.1e2 

2.0e2 

1.8e3 

1.2e2 

7.0e2 

l.le2 

4.4el 

1 

3 

1.8e-3 

2.5e3 

3.6e3 

3.0e4 

4.7e2 

9.4e3 

3.1e3 
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1 

4 
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3 
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1 
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3 
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3.1e-6 
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7.8e4 

5.3e3 

2.5e3 
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6 

10 

4.3e-7 

2.0e8 

3.1e8 

2.0e9 

1.7e6 

9.6e4 

2.3e4 

1.6e4 

6 

11 

4.3e-9 

2.el0 

4.el0 

4.ell 

5.7e8 

3.7e5 

3.9e5 

2.6e5 
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Table  2:  Condition  numbers  for  problem  afiro 
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Figure  2:  Eigenvalues  for  A'2,  A'3,  A’4  for  problem  a/iro 
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7.  Conclusions 

For  symmetric  indefinite  systems  of  linear  equations,  we  have  shown  that  the  Bunch- 
Parlett  factorization  can  be  used  to  provide  a  preconditioner  for  the  Paige-Saunders 
algorithm  SYMMLQ  (Section  2).  This  general  result  led  us  to  consider  iterative 
methods  for  the  KKT  systems  arising  in  barrier  methods  for  QP  and  nonlinear 
programming.  The  preconditioner  (3.7)  should  play  an  important  role  in  future 
interior- point  implementations  for  large-scale  constrained  optimization. 

For  linear  programs,  the  sensitivity  analysis  associated  with  the  partitioned  KKT 
system  (3.5)  led  us  to  consider  the  true  sensitivity  of  A',  as  reflected  by  the  pre¬ 
conditioner  and  the  transformed  system  Kx  (4.1)-(4.2).  In  turn,  the  fact  that 
cond(A'j)  is  typically  much  smaller  than  cond(AD2AT)  motivated  development  of 
the  preconditioners  M2,  M3,  M4  (4.3),  (4.5),  (4.7). 

Subject  to  effective  methods  for  choosing  B,  we  expect  these  KKT  precondition¬ 
ers  to  bring  improved  reliability  to  interior-point  LP  algorithms.  Implementations 
based  on  direct  or  iterative  solves  with  AD2AT are  often  remarkably  effective,  but  the 
extreme  ill-conditioning  of  the  AD2AT  systems  as  solutions  are  approached  makes 
their  use  tantamount  to  walking  the  razor’s  edge. 

A  switch  to  the  full  KKT  system  should  be  beneficial  as  Gay  [Gay89]  sug¬ 
gests,  particularly  when  A  contains  some  relatively  dense  columns  that  prevent 
exact  Cholesky  factorization  of  AD2AT.  Fortunately,  since  B  becomes  more  sharply 
defined  near  a  solution,  the  KKT  preconditioners  will  become  most  effective  when 
they  are  most  needed. 
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A.  A  preconditioner  from  the  Bunch-Parlett  factorization 

The  following  Fortran  77  routine  illustrates  the  construction  of  a  positive-definite 
matrix  M  =  UTDU  from  the  Bunch-Parlett  factorization  A  =  UTDU  produced  by 
the  Harwell  MA27  package  of  Duff  and  Reid  [DR82,DR83]. 

Subroutine  syprec  overwrites  the  representation  of  D  in  the  MA27  data  struc¬ 
ture.  A  typical  application  would  contain  calls  of  the  form 

call  ma27ad(  n,  nz,  ...  ) 
call  ma27bd(  n,  nz,  ...  ) 
call  syprecC  n,  la,  ...  ) 

to  factorize  A  and  compute  D,  followed  by  multiple  calls  of  the  form 
call  ma27cd(  n,  a  ,  ...  ) 
to  solve  systems  involving  M. 

*  - 

subroutine  syprecC  n,  la,  liw,  a,  iw,  negl,  neg2  ) 

implicit  double  precision  (  a-h,  o-z  ) 

double  precision  a(la) 

integer* 2  iw(liw) 

integer  negl,  neg2 


syprec  (symmetric  PREConditioner)  takes  the  factors 
A  =  U*  D  U 

from  Duff  and  Reid's  Harwell  subroutine  KA27BD  and  changes  the 
block-diagonal  matrix  D  to  be  a  positive-definite  matrix  Dbar  with 
the  same  lxl  and  2x2  block-diagonal  structure. 

The  eigensystem  D  =  Q  E  Q’  is  used  to  define  Dbar  =  Q  lEl  Q’, 
where  I  El  contains  the  absolute  values  of  the  eigenvalues  of  D. 

The  matrix 

Abar  =  U'  Dbar  U 

is  then  an  exact  preconditioner  for  A,  in  the  sense  that  SYMMLQ 
would  take  only  2  iterations  to  solve  Ax  =  b  Cor  1  iteration  if 
D  =  Dbar  is  already  positive  definite). 

If  the  original  matrix  A  is  close  to  some  other  matrix  K, 

Abar  should  be  a  good  preconditioner  for  solving  K  x  *  b. 

Mote  that  MA27  stores  the  elements  of  D(inverse)  and  (  -  U  ) 
within  A  and  IV.  However,  modifying  a  2x2  block  of  D(inverse) 
looks  the  same  as  modifying  the  2x2  block  itself. 


10  Mar  1989:  First  version. 

Systems  Optimization  Laboratory,  Stanford  University. 
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intrinsic  abs  ,  sqrt 

integer  alen,  apos 

logical  single 

parameter  (  zero  =  0.0d+0,  one  =  1.0d+0,  two  =  2.0d+0  ) 

negl  =  0 

neg2  =  0 

nblk  =  abs(  iw(l)  ) 

ipos  =  2 

apos  =  1 

do  100,  iblk  =  1,  nblk 
ncols  -  iw(ipos) 

if  (ncols  .It.  0)  then 
nrows  =  1 

ncols  =  -  ncols 
else 

ipos  »  ipos  +  1 
nrows  =  in (ipos) 
end  if 

*  Process  the  diagonals  in  this  block. 

alen  ■  ncols 
single  =  .true. 

do  SO,  k  =  ipos  +  1,  ipos  +  nrows 
if  (  single  )  then 
alpha  =  a(apos) 

j  =  iw(k) 

single  =  j  ,gt.  0 

if  (  single  )  then 

if  (  alpha  .It.  zero  )  then 


The  lxl  diagonal  is  negative. 


negl  =  negl  +  1 
a(apos)  =  -  alpha 
end  if 

else 

beta  ■  a(apos+l  ) 
gamma  =  a(apos+alen) 

if  (  alpha  *  gasuaa  .It.  beta**2  )  then 


The  2x2  diagonal  is  indefinite. 
Find  its  eigensystea  in  the  fora 


* 

* 

e 

* 


* 

* 

* 
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(  alpha  beta  )  =  (  c  s  )  (  «1  )  (  c  s  ) 

(  beta  gamma  )  (s-c)(  e2)(s-c) 


tau 

= 

(  gamma  -  alpha  )  /  (  two  *  beta  ) 

t 

= 

abs(  tau  )  +  sqrt(  tau*<*2  +  one  ) 

t 

= 

-  one  /  t 

if  ( 

tau 

.It.  zero  )  t  =  -  t 

c 

= 

one  /  sqrtC  t**2  +  one  ) 

s 

= 

t  *  c 

el 

- 

alpha  *  beta  *  t 

e2 

= 

gamma  -  beta  *  t 

Change  el  and  e2  to  their  absolute  values 
and  then  multiply  the  three  2x2  matrices 
to  get  the  modified  alpha,  beta  and  gamma. 


if  (  el  .It.  zero  )  then 
neg2  =  neg2  +  1 

el  =  -  el 

end  if 

if  (  e2  .It.  zero  )  then 
neg2  ■  neg2  +  l 


e2  =  -  e2 
end  if 

alpha  = 

beta  = 

gamma  = 

aCapos  )  = 

a(apos+l  )  = 

a(apos+alen)  = 

end  if 
end  if 
else 

single  =  .true, 
end  if 


c**2 

*  el 

+  s**2  * 

e2 

c*s 

*(el 

-  e2) 

s**2 

*  el 

♦  c**2  * 

62 

alpha 

beta 

gamma 


apos  =  apos  +  alen 

alen  =  alen  -  1 

50  continue 


ipos  =  ipos  +  ncols  +  1 
100  continue 


*  end  of  syprec 
end 
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We  discuss  the  solution  of  sparse  linear  equations  Ky  —  z,  where  K  is  symmetric  and  indefinite. 
Since  exact  solutions  are  not  always  required,  direct  and  iterative  methods  are  both  of  interest. 

An  important  direct  method  is  the  Bunch-Parlett  factorization  K  —  UT DU ,  where  U  is  triangular 
and  D  is  block-diagonal.  A  sparse  implementation  exists  in  the  form  of  the  Harwell  code  MA27  An 
appropriate  iterative  method  is  the  conjugate-gradient-like  algorithm  SYMMLQ,  which  solves  indefinte 
systems  with  the  aid  of  a  positive-definite  preconditioner. 

For  any  indefinite  matrix  K ,  we  show  that  the  Ur DU  factorization  can  be  modified  at  nominal 
cost  to  provide  an  “exact”  preconditioner  for  SYMMLQ.  We  give  code  for  overwriting  the  block-diagonal 
matrix  D  produced  by  MA27. 

We  then  study  the  KKT  systems  arising  in  barrier  methods  for  linear  and  nonlinear  programming, 
and  derive  preconditioners  for  use  with  SYMMLQ. 

For  nonlinear  programs  we  suggest  a  preconditioner  based  on  the  “smaller"  KKT  system  associated 
with  variables  that  are  not  near  a  bound.  For  linear  programs  we  propose  several  preconditioners  based 
on  a  square  nonsingular  matrix  B  that  is  analogous  to  the  basis  matrix  in  the  simplex  method  The  aim 
is  to  facilitate  solution  of  full  KKT  systems  rather  than  equations  of  the  form  AD2  An  =  r  when  the 
latter  become  excessively  ill-conditioned. 
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